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Abstract 

The semi-classical Wigner-Kirkwood h expansion method is used to calculate shell corrections for 
spherical and deformed nuclei. The expansion is carried out up to fourth order in h. A systematic 
study of Wigner-Kirkwood averaged energies is presented as a function of the deformation degrees 
of freedom. The shell corrections, along with the pairing energies obtained by using the Lipkin- 
Nogami scheme, are used in the microscopic-macroscopic approach to calculate binding energies. 
The macroscopic part is obtained from a liquid drop formula with six adjustable parameters. 
Considering a set of 367 spherical nuclei, the liquid drop parameters are adjusted to reproduce 
the experimental binding energies, which yields a rms deviation of 630 keV. It is shown that the 
proposed approach is indeed promising for the prediction of nuclear masses. 
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I. INTRODUCTION 



Production and study of loosely bound exotic nuclei using Radioactive Ion Beam facilities 
is of current interest Q, O] . These experiments have given rise to a number of interesting and 
important discoveries in nuclear physics, like neutron and proton halos, thick skins, disap- 
pearance of magicity at the conventional numbers and appearance of new magic numbers, 
etc. Further, advances in detector systems, and in particular, the development of radioactive 
beam facilities like Spiral, REX-Isolde, FAIR, and the future FRIB may allow to investigate 
new features of atomic nuclei in a novel manner. 

The study of nuclear masses and the systematics thereof is of immense importance in 
nuclear physics. With the advent of mass spectrometry, it is possible to measure masses of 
some of the short lived nuclei spanning almost the entire periodic table 0, 4]. For example, 
;he ISOL (isotope separator online) based mass analyzer for superheavy atoms (MASHA) 
5), [fj] coming up at JINR-Dubna will be able to directly measure the masses of separated 
atoms in the range 112 < Z < 120. The limitation on measurements is set by the shortest 

n 2 

measurable half-life, Ti/ 2 ~ 1.0 s [5[. The JYFLTRAP 7] developed at the University of 
Jyvaskyla, on the other hand, enables to measure masses of stable as well as highly neutron 
deficient nuclei (for masses up to A = 120) with very high precision (~50 keV) 7|. 

On the theoretical front as well, considerable progress has already been achieved in the 
accurate prediction of the nuclear masses, and it is still being pursued vigorously by a number 
of groups around the globe. This is of great importance, since an accurate knowledge of the 
nuclear masses plays a decisive role in a reliable description of processes like the astrophysical 



r-process (see, for example, 



There are primarily two distinct approaches to calculate 



masses: a) the microscopic nuclear models based on density functional theory like, Skyrme 
8, 9] and Gogny [k]] Hartree-Fock-Bogoliubov or Relativistic Mean Field (RMF) models 
ll|), b) microscopic-macroscopic (Mic-Mac) models 



The Mic-Mac models are based on the well-known Strutinsky theorem. According to this, 
the nuclear binding energy, hence the mass can be written as sum of a smooth part, and 
an oscillatory part which has its origins in the quantum mechanical shell effects. The latter 
consists of the shell correction energy and the pairing correlation energy which in the Mic- 
Mac models are evaluated in an external potential well. The smooth part is normally taken 
from the liquid drop models of different degrees of sophistication. The largest uncertainties 
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arise in the calculation of shell corrections. The shell correction is calculated by taking 
the difference between the total quantum mechanical energy of the given nucleus, and the 
corresponding 1 averaged ' energy. Usually, the averaging is achieved by the well-established 
Strutinsky scheme 



[la. ~v\ . This technique of calculating the averaged energies runs into 
practical difficulties for finite potentials, since for carrying out the Strutinsky averaging, one 
requires the discrete single-particle spectrum, with cut-off well above (at least 3ftwo, 
being the major shell spacing) the Fermi energy. For a realistic potential, this condition is 
not met, since continuum may start within ~ hwo of the Fermi energy. Standard practice is 
to discretise the continuum by diagonalising the Hamiltonian in a basis of optimum size. A 



number of Mic-Mac calcu 
(see, for example 
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ations with varying degree of success are available in the literature 



Il5l|). The Mic-Mac models typically yield better than ~0.7 
MeV rms deviation in the masses. All these models agree reasonably well with each other 
and with experiment, but deviate widely among themselves in the regions far away from the 
valley of stability. 
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25|, on 



The semi-classical Wigner-Kirkwood (WK) approach [18, 
the other hand, makes no explicit reference to the single-particle spectrum, and achieves an 
accurate averaging of the given one-body Hamiltonian. Thus, the WK approach is a good 
alternative to the conventional Strutinsky smoothing scheme. The quantum mechanical 
energy is calculated by diagonalising the one-body Hamiltonian in the axially symmetric 
deformed harmonic oscillator basis with 14 shells. The difference between the total quantum 
mechanical energy and the WK energy in the external potential well yields the value of the 
shell correction for a given system. In the present work, we propose to carry out a reliable 
microscopic-macroscopic calculation of the nuclear binding energies (and hence the masses) 



employing the semi-classical Wigner-Kirkwood (WK) h expansion [18 



m 
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25| for the calculation of shell corrections instead of the Strutinsky scheme. An exploratory 



study of using the WK method to compute the smooth part of the energy has been reported 
earlier to test the validity of the Strutinsky scheme, especially near the driplines 271 ] . 

It is known that the WK level density {gwK{s)) with the h 2 correction term exhibits a 
£-V 2 divergence as e — » 0, for potentials which vanish at large distances as for instance 
Woods-Saxon potentials (see, for example, Ref. |26|). The Strutinsky level density, on the 



28] 



contrary, exhibits only a prominent peak as e — > 0. It was therefore concluded in Ref. 
that the divergence of the WK level density as e — > is unphysical, and the Strutinsky 
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smoothed level density should be preferred. It should however be noted that the WK 
level densities, energy densities, etc., have to be understood in the mathematical sense of 
distributions and, consequently, only integrated quantities are meaningful. In fact, it has 
been shown |25| that the integrated quantities such as the accumulated level densities are 
perfectly well behaved, even for e — > 0. 

Pairing correlations are important for open shell nuclei. In the present work, these are 
;aken into account in the approximate particle number projected Lipkin-Nogami scheme 



29l . l30l . l3l| . Odd-even and odd-odd nuclei are treated in an entirely microscopic fashion 
(odd nucleon blocking method in the uniform filling approximation), allowing an improved 
determination of odd-even mass differences, see e.g. the discussion in [321 ] . The majority 
of nuclei in the nuclear chart are deformed. In particular, it is well known that inclusion 
of deformation is important for reliable predictions of nuclear masses. Therefore, here we 
incorporate in all three deformation degrees of freedom (/3 2 , At, 7). To our knowledge, 
no such detailed and extensive calculation based on the WK method is available in the 
literature. 

The paper is organised as follows. We review the WK expansion in Section 2. The choice 
of the nuclear, spin-orbit, and Coulomb potentials forms the subject matter of Section 3. 
Details of the WK calculations are discussed in Section 4. A systematic study of the WK 
energies for neutrons and protons as a function of the deformation degrees of freedom is 
presented in Section 5. The shell corrections for the chains of Gd, Dy and Pb isotopes 
obtained by using our formalism are reported, and are compared with those calculated 
employing the traditional Strutinsky averaging technique, in Section 6. Section 7 contains 
a brief discussion on the Lipkin-Nogami pairing scheme. As an illustrative example, the 
calculation of the binding energies for selected 367 spherical nuclei is presented and discussed 
in Section 8. Section 9 contains our summary and future outlook. Supplementary material 
can be found in appendices A and B. 

II. SEMI-CLASSICAL WIGNER-KIRKWOOD EXPANSION 



Following Ref . 20J| , we consider a system of N non- interacting fermions at zero temper- 
ature. Suppose that these fermions are moving in a given one-body potential including the 
spin-orbit interaction. To determine the smooth part of the energy of such a system, we 
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start with the quantal partition function for the system: 



Z(/3) = Tr(exp (-/?£)). 



Here, H is the Hamiltonian of the system, given by: 

-h 2 



H 



2m 



V 2 + V(f) + V LS {r) 



(1) 



(2) 



where V{r) is the one-body central potential and Vls(t) ls the spin-orbit interaction. 

In order to average out shell effects, the simplest one could do is replace the partition 
function in the above expression by the classical partition function. That is, one replaces 
the Hamiltonian in Eq. (OQ) by the corresponding classical Hamiltonian. This yields the 
well-known Thomas- Fermi equations for particle number and energy. Way back in 1930 's, 
E. Wigner [l8| and J. G. Kirkwood [3] developed a systematic expansion of the partition 
function in powers of the Planck's constant, h, its first term bein g th e classical partition 



function. Details of this method can be found in Refs. l2lL I22L |23| . Such expansion 

of the quantal partition function in powers of H is often known as Wigner-Kirkwood (WK) 
expansion. Systematic corrections to the Thomas-Fermi energy can be obtained by using 
the WK expansion. 

In this work, we shall use the WK expansion up to fourth order. For brevity, we represent 
the potentials and form factors without mentioning the dependence on the position vector. 
Ignoring the spin-orbit interaction, the WK expansion of the partition function, correct up 
to fourth order is given by [20] ]: 

(3 2 h 2 



4tt 3 / 2 \h 2 J 

33 / j-9 \ 2 



1 - 



24m 



-V Z V 



+ 



A 



1440 



(2^) {- 7v4y + 5 ^( vV ) 2 + Pv 2 (vvy 



(3) 



The spin-orbit interaction, in general, can be written as: 

lkH 2 



V, 



LS 



2m 



V/xVU 



(4) 



where a is the unit Pauli matrix, k is the strength of spin-orbit interaction, and / is the 
spin-orbit form factor. With the inclusion of such spin-orbit interaction, the WK expansion 
for the full partition function splits up into two parts: 



7 (4) 
WK 



((3) = + Z^((3) . 



(5) 



Here, Z^(P) is given by Eq. ([3]), and the spin-orbit contribution to the partition function, 
correct up to fourth order in H, reads [20I ]: 

k 2 ^/ 2 [2m\ 1/2 



where 



7 (4) 

LS 87T 3 / 2 



+ 



/3 1 / 2 /fr 2 \ 



96vr 3 / 2 



h = -P (v/) 2 (W) + ^V 2 (V/) 2 - (v 2 /) 2 + V/ • V (V 2 /) 



h = (v/rv 2 /-2v/-v(v/r 

h = (V/) 4 . 



(6) 



(7) 

(8) 
(9) 



The level density gwK, particle number N, and energy E can be calculated directly from 
the WK partition function by Laplace inversion: 



g WK (e) = £?Z$ K {0) 



(10) 



N = C 



-1 { Zwk(P) 
P 



(11) 



and 



-1 ( Zwk(P) 

a 2 



E 



\N - C 



(12) 



where A is the chemical potential, fixed by demanding the right particle number, and 
denotes the Laplace inversion. Using the identity 



(a - vy - 1 



-S(X-V) ,for/i > 



and noting that, in order to get inverse Laplace transforms in convergent form, 

-1 de'W 



-pv 



(3 dV ' 

one obtains the level density for each kind of nucleons assuming spin degeneracy: 



9WK\£) 



1 / if 



3tt 2 \h 2 



l (f - l ' )1/2 + ^{i K2(v/ ) 2 ( e -^ 1/2 



+-AV(e-Vy 3 ' 2 
lb 



(13) 



(14) 



e(e-V), (15) 



the particle number: 
N 



1 < 2m \ m I ir 



3tt 2 \h 2 J 



and the energy: 

E = XN - XC 2 ^' 2 Idr 



(a - vo 3/2 - — (a - vr 1/2 w 

32m 

+ 3 -^(X-Vf 2 (Vf) 2 



6m 



e(\-v), (16) 



3tt 2 \h 2 J 

i /^y /2 



|(A-y) 5/2 --^(A-y) 1/2 W 
5 lorn 



e (A - v) 



5760tt 2 V2m/ 



J df(\-vy 1/2 {-7vV} 

_i y df (A - v)- 3 / 2 {5 (vV) 2 + V 2 (W) 2 } 



e (A - v) 



67T 2 



^ Jdr(X-Vf 2 (Vffe(X-V) 



1 (— I df(X-V) 1/2 



^ 2 {^v 2 (v/) 2 -(v 2 /) 2 + v/-v(v 2 /) 



48tt 2 V2my 

~ ^A^f [ " 2 ^ i (V/r VV " 2 V/ ' V (V/r ^ + 2/?4 (V/r 



e (a - v) 

(17) 



It should be noted that we have explicitly assumed that all the derivatives of the potential 
V and the spin-orbit form factor / exist. The expansion defined here is therefore not valid 
for potentials with sharp surfaces. This automatically puts a restriction on the choice of 
the Coulomb potential: the conventional uniform distribution approximation for the charge 
distribution cannot be used in the present case. We shall discuss this point at a greater 
length in the next section. The integrals in the above expressions are cut off at the turning 
points, defined via the step function. The chemical potential A appearing in these equations 
is determined from Eq. (ITBT) . separately for neutrons and protons. Further, it is interesting 
to note that the spin-orbit contribution to the particle number iV as well as to the energy E 
appears only in the second order in h. Secondly, the level density and particle number are 
calculated only up to the order h 2 . It can be shown 20j that for the expansion correct up to 
fourth order in h, it is sufficient to take Z^) K up to order h 2 in Eq. (fTTT) to find the chemical 
potential (and hence the particle number), whereas one has to take the full partition function 
Zyy K up to order h 4 in Eq. (T121) to compute the energy in the WK approach. 
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The divergent terms appearing in Eq. §1 
the chemical potential. Explicitly: 

d x (A - V) 1/2 
d\ (A - Vf 2 



are treated by differentiation with respect to 



i ( A-vr 1/2 



-i(X-V) 



-3/2 



(18) 
(19) 



In practice, the differentiation with respect to chemical potential is carried out after evalua- 
tion of the relevant integrals. Numerically, this approach is found to be stable. Its reliability 
has been checked explicitly by reproducing the values of fourth-order WK corrections quoted 
in Ref. Q- 

The WK expansion thus defined, converges very rapidly for the harmonic oscillator po- 
tential: the second-order expansion itself is enough for most practical purposes. The conver- 
gence for the Woods-Saxon potential, however, is slower than that for the harmonic oscillator 
potential, but it is adequate 33j. For example, for ~ 126 particles, the Thomas- Fermi en- 
ergy is typically of the order of 10 3 MeV, the second-order (h 2 ) correction contributes a few 
10's of MeVs, and the fourth-order (^i 4 ) correction yields a contribution of the order of 1 
MeV. This point will be discussed in greater details later. It is also important to note that 
the WK h expansion of the density matrix has a variational character and that a variational 
theory based on a strict expansion of the of h has been established 341 ]. 

The WK approach presented here should be distinguished from the extended Thomas- 
Fermi (ETF) approach. Divergence problems at the classical turning points (see the particle 
number and energy expressions above) can be eliminated by expressing the kinetic energy 
density as a functional of the local density. This is achieved by eliminating the chemical 
potential, the local potential, and the derivatives of the local potential (for further details, 



see Ref. 



351]). It cannot be accomplished in closed form, and has to be done iteratively, 



leading to a functional series for the kinetic energy density. The resulting model is what 
is often referred to as the ETF approach. The WK approach as presented here, in this 
sense, js^ the starting point for ETF approach (further details of ETF can be found in 



Refs. 



22 



23 



25 



36 



33, 



38j). The conventional ETF approach exhibits somewhat slower 



convergence properties w' 
each given power in h 125 



rich has been attributed to a non-optimal sorting out of terms of 



35]. 
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III. CHOICE OF POTENTIAL 



Form of the Nuclear Potential 



The spherically symmetric nuclear mean field is well represented by the Woods-Saxon 
(WS) form [39[], given by: 

Vo 



V(r) 



(20) 



1 + exp ((r - Ro)/a) 

where Vo is the strength of the potential, Rq is the half-density radius, and a is the dif- 
fuseness parameter. The WS form factor defined here, can be easily generalised to take the 
deformation effects into account. Note that the distance function l(r) = r — R appearing in 
Eq. (1201) can be interpreted as the minimum distance of a given point to the nuclear surface, 
defined by r = R . One might thus generalise it to the case of deformed surfaces as well. 
Using the standard expansion in terms of spherical harmonics, a general deformed surface 
may be defined by the relation r = r s , where 



CR (1 + J>V y A, A 

A, /Li 



(21) 



Here, the functions are the usual spherical harmonics and the constant C is the volume 
conservation factor (the volume enclosed by the deformed surface should be equal to the 
volume enclosed by an equivalent spherical surface of radius R ): 

-1/3 

(22) 



C 



The distance function to be used in the WS potential would be the minimum distance of 
a given point to the nuclear surface defined by r = r s . Such definition has been used quite 



extensively in the literature, with good success (see, for example, Refs. (40 



41 



3,Q, 



44j) 



However, in the present case, this definition is not convenient, since it should be noted that 
the calculation of this distance function involves the minimisation of a segment from the 
given point to the nuclear surface. This in turn implies that each calculation of the distance 
function (for given r, 9, and <p coordinates: we are assuming a spherical polar coordinate 
system here) involves the calculation of two surface angles 9 S and <f) s , and these are implicit 
functions of r, 9, and 0. See Fig. (8) in Appendix A for details. Since the WK calculations 
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involve differentiation of the WS function, one also needs to differentiate 9 S and <f> a , which 
are implicit functions of r, 8, and </>. 

Alternatively, the distance function for the deformed Woods-Saxon potential can be writ- 
ten down by demanding that the rate of change of the potential calculated normal to the 
nuclear surface and evaluated at the nuclear surface should be a constant [45| which, indeed, 
is the case for the spherical Woods-Saxon form factor. Thus, 

{h- W(f)} r=rs = constant, (23) 

where n is the unit vector normal to the surface (r = r s ) and is given by 

_ V(r - r a ) 

~ |V(r - r s )\ ■ {2A} 
In fact, the above condition (T2"3"j) is related to the observation that the second derivative of 
the spherical Woods-Saxon form factor vanishes at the nuclear surface, defined by r = Rq. 
The resulting distance function is given by 46J: 



= M^fc ' (25) 

where r s is as defined in Eq. ( l2Tj) . The denominator is evaluated at r = r s . Writing the 9 
and <p derivatives of r s as A and B respectively, we get: 

= (r ~ rs) , (26) 
A/1 + l 2 /r s 2 

with 

7 2 = A 2 + EWtf . (27) 

In the present work, we use the distance function as defined in Eq. ( T25l) . The WS potential 
thus reads 

V (r) = Yl . (28) 

{) l + exp(/(fO/a) { ! 

It is straightforward to check that the Woods-Saxon potential defined with the distance 
function as given by Eq. (f2"5|) satisfies the condition ([231 . Substituting this Woods-Saxon 
potential in h ■ W (r), we get 

h ■ W (r) = —f(r) (f(r) - l)n-Vl(f) 
a 

= -f(?) (/(rl - 1) 
a 

+ (r - r - ) W^T- V |V(r-r,)U. ]-( 2i " 
10 





V (r - r s ) 




V (r - t 8 ) 


r- 





Here, f(r) = [1 + exp (l{f)/a)\ 1 is the Woods-Saxon form factor. Clearly, at the surface 
defined by r = r s , the quantity h ■ W (f) is constant. 



B. Deformation Parameters 



In practice, we consider three deformation degrees of freedom, namely, fa, fa and 7. 
These parameters are related with the parameters a\ jfl introduced in Eq. (l2~Tj) . Note that 
for the given choice of deformation parameters, A takes values 2 and 4. The projection \i 
takes the values 0, ±2 for A = 2 and the values 0, ±2 and ±4 for A = 4. Further, existence 



of symmetry planes (x,y), (y,z) and (z,x) implies that [42] 



«2,2 = «2 -2, «4,2 = a4,-2, «4,4 = «4,-4 • 

Thus, we get: 

r s (6, <P) = CR [1 + « 2i0 F 2 ,o (6) + a 2 , 2 {Y 2 , 2 (6, 0) + F 2 ,- 2 (0, 0)} + 04,0*4,0 (6) 

+ a 4 ,2 {^4,2 (0, <P) + ^4,-2 (0, <P)} + 04,4 {^4,4 (0, 0) + ^4,-4 (0, <P)}] , 

where 



02,0 

«2,2 
"4,0 
«4,2 
O4 4 



fa cos 7 



6 



2^2 sin 7 
/3 4 (5 cos 2 7 + 1) 



-\/^4sm2 7 



For further details, see Ref. 



42]. 



70 



144 



fa sin 2 7 . 



(30) 

(31) 
(32) 
(33) 

(34) 

(35) 



C. Woods- Saxon Parameters 



The parameters 47] appearing in the Woods-Saxon potential are as defined below: 

1. Central potential: 
a. Strength: 

N-Z) , x 

Ui = -'oilT^— (36) 
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with [7o=53.754 MeV and f7 a =0.791. 
b. Half-density radius: 



N — Z ' 

R = roA 1 / 3 { 1 T Ci— j— \ + c 2 (37) 



.4 

with r =1.19 fm, Ci=0.116 and C2=0.235 fm. 

c. Diffuseness parameter: assumed to be same for neutrons and protons, and has the 
value a = 0.637 fm. 

2. Spin-orbit potential: 

a. Strength: 

ft 2 ( N — 7\ 

Vso = X U — 2 {l*U 2 — j (38) 

with f7o=53.754 MeV, A =29.494 and C/ 2 =0.162. 

b. Half density radius and diffuseness parameter are taken to be the same as those for 
the central potential. 

In these expressions, the + (— ) sign holds for protons (neutrons). 

The parameters have the isospin dependence of the central and spin-orbit potentials 
"built-in". This potential yields a reasonably good description of charge radii (both mag- 
nitude and isospin dependence) as well as of moments of inertia for a wide range of nuclei. 
It has been used extensively in the total Routhian surface (TRS) calculations, and it has 
been quite successful in accurately reproducing energies of single-particle as well as collective 



states 



48|]. 



D. Coulomb potential 

The Coulomb potential is calculated by folding the point proton density distribution 
p(r*), assumed to be of Woods-Saxon form. For simplicity, its parameters are assumed to be 
the same as those for the nuclear potential of protons. The reason for using folded potential 
here is, as we have indicated in section II, the WK expansion is not valid for potentials with 
sharp surfaces. 

The Coulomb potential for the extended charge distribution is given by: 

V c (r) = e 2 J p(f)^—df . (39) 

12 



Here, 



\f-r\ = {r 2 + r' 2 - 2rr'cos\^} 1/2 , (40) 

where 

cos^ = cos ^ cos 6 1 ' + sin # sin cos(0 — <fi') , (41) 
as explained in Appendix A. 
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FIG. 1: Coulomb potentials obtained by using diffuse density and sharp surface approximation for 
208 Pb. 



It is instructive at this point, to compare the Coulomb potential calculated from the 
diffuse density with the corresponding potential obtained by using the conventional uniform 
density (sharp surface) approximation. Such comparison for 208 Pb is plotted in Fig. [TJ The 
radius parameter for the diffuse density approach as well as for the sharp surface approxima- 
tion is assumed to be equal to 7.11 fm (see the discussion on the choice of the Woods-Saxon 
parameters in Section 3). It can be seen that in the exterior region, the two potentials agree 
almost exactly, as expected. In the interior, however, the potential obtained from the dif- 
fuse density turns out to be somewhat less repulsive than that from the density with sharp 
surface. 
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IV. DETAILS OF THE WK CALCULATIONS 



In the present work, we restrict our calculations to three deformation degrees of freedom, 
namely, fa, fa and the angle 7. The inclusion of 7 allows to incorporate triaxiality. Thus, 
the present WK calculation is genuinely three dimensional. In principle it is natural to use 
a cylindrical coordinate system here. The spherical polar coordinates, however, turn out to 
be more convenient. The reason is, the cylindrical coordinates involve two length variables, 
and one angular coordinate which means that the turning points have to be evaluated for 
two coordinates (p and z). This makes the calculations very complicated. On the other 
hand, the spherical polar coordinates involve only one length variable, and thus the turning 
points are to be evaluated only for one coordinate (r). The numerical integrals involved are 
evaluated using Gaussian quadrature. 

The first step in the WK calculations is the determination of the chemical potential. 
This has to be done iteratively, using Eq. ffl6|) . Since the turning points are determined 
by the chemical potential, they have to be calculated using a suitable numerical technique 
at each step. Once the values of the chemical potential are known, the WK energies up 
to second order can be calculated in a straightforward way. The fourth-order calculations 
are very complicated, since they require higher-order derivatives of nuclear potentials, spin- 
orbit form factors, and the Coulomb potential. The former can be evaluated analytically 
in the present case. The expressions are extremely lengthy, and we do not present them 
here. Comparatively, the derivatives of the Coulomb potential look simple; the Laplacian 
and Laplacian of Laplacian are completely straightforward: the former is proportional to 
the proton density and the latter is just the Laplacian of the WS form factor. However, the 
calculations also need terms like Laplacian of the gradient squared of the total potential. In 
the case of protons, this involves one crossed term: 



where Vc is the Coulomb potential and Vn is the nuclear potential. The determination of 
such objects is tricky. It turns out that if one uses the form of the Coulomb potential defined 
above, the calculation of expression (142 p becomes numerically unstable. 
There exists an alternative for of the Coulomb potential: 



V 2 (W c (r) ■ VVy(rO) 



(42) 




(43) 
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where the notation Vj» means that the Laplacian is calculated with respect to the variables 
r', 9', and (ft'. Eqs. (139]) and fH3|) are exactly equivalent. This is proved explicitly in 
Appendix B. With this form, one can calculate the first and second derivatives (not the 
Laplacian) of the Coulomb potential. Calculation of the higher-order derivatives of the 
Coulomb potential, even with the form defined in Eq. (1431) . turns out to be numerically 
unstable. For this purpose, we employ the Poisson's equation. The details are presented in 
Appendix B. Once all the derivatives are known, the fourth-order WK calculations can be 
carried out. 

It turns out that the WK calculations for protons are very time consuming. This is 
due to the fact that the calculation of Coulomb potential (Eq. (I39p ). in general, involves 
evaluation of three dimensional integral for each point (r,6,<f)). Typically, it takes few tens 
of minutes to complete one such calculation. This is certainly not desirable, since our aim 
is to calculate the masses of the nuclei spanning the entire periodic table. To speed up the 
calculations, we use the well-known technique of interpolation. Since we are using spherical 
polar coordinates, the turning points are to be evaluated only for the radial coordinate, r. 
For the entire WK calculation, the 9 and mesh points remain the same (over the domains 
[0, 7r] and [0, 2tt], respectively), whereas the r mesh points change from step to step. This 
happens in particular during the evaluation of the chemical potential. Once the convergence 
of the particle number equation (Eq. ffTB"]) ) is achieved, the r mesh points as well, remain 
fixed. 

Motivated by the above observations, we apply the following procedure: 

1. Before entering into the actual WK calculations (determination of chemical potential, 
etc.), for each pair of 9 and mesh points, we calculate the Coulomb potential (Eq. 
(1391) ) over a range to 16 fm at equidistant radial mesh points (the typical mesh size 
being 0.1 fm). 

2. Next, for each pair of 9 and mesh points, we fit a polynomial of degree 9 in the 
radial coordinate r to the Coulomb potential calculated in the above step. Thus, the 
fitting procedure is to be repeated N e x times, N e (N^) being total number of mesh 
points for the 9 ((f)) integration. 

Thus, for any given value of radial the coordinate r (and fixed 9 and 0), the Coulomb 
potential can be easily calculated just by evaluating the 9 th degree polynomial in r. It is 
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found that this interpolation procedure is very accurate. The maximum percentage difference 
between the fitted and the exact Coulomb potentials is 0.4% for a highly deformed nucleus. 

V. VARIATION OF WIGNER-KIRKWOOD ENERGIES WITH DEFORMATION 
PARAMETERS 

A sample WK calculation is performed for system of 126 neutrons and 82 protons. The 
variation of the Thomas-Fermi energy and of the different correction terms as a function of 
the quadrupole deformation parameter 02 is presented in Fig. (j2J). The other two defor- 
mation parameters, and 7, are set to zero in this test case. The partial contributions 
to the WK energy are plotted separately for protons and neutrons. It is found that all the 
correction terms vary smoothly as a function of deformation. As expected, the value of the 
contributions from the h 2 and h 4 terms to the averaged energy decreases rapidly. It is found 
that the proton and neutron Thomas-Fermi energies have opposite trends with respect to 
increasing j3 2 - If Coulomb potential is suppressed, it is found that the Thomas- Fermi ener- 
gies for protons follow the same trend as those for the neutrons. Further, it is interesting 
to note that comparatively, the variation in the second-order corrections with respect to 
deformation parameters is stronger than that in the Thomas-Fermi energies (~ 10% for 
second-order corrections and ~ 3% for Thomas- Fermi energies). 

Next, the variation of the Thomas- Fermi energy and of the correction terms as a function 
of the hexadecapole deformation parameter /3 4 is plotted in Fig. Q. Here, /3 2 is taken to 
be 0.2 and 7 is set to zero. It is seen that again, the different energies vary smoothly as a 
function of (3^. The Thomas- Fermi energy for protons is found to have very little variation 
with respect to the (3^ deformation parameter. In contrast, the corresponding energies 
for neutrons have a stronger dependence on /5 4 . The same behaviour is also observed in 
the corresponding quantum mechanical energies. It is found that the proton and neutron 
Thomas- Fermi energies have a very similar behaviour if the Coulomb potential is suppressed. 
Further, to check if this conclusion depends on the value of j3 2 , the analysis is repeated for 
P2 = 0.4, and the same conclusion is found to emerge. 

The behaviour of the Thomas-Fermi energies for protons in the above cases (Figs. (2) 
and (3)) seems to be due to the Coulomb potential. In the case of variation with respect to 
/3 2 , qualitatively it can be expected that with increasing quadrupole deformation, protons 
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FIG. 2: Wigner Kirkwood energies for 126 neutrons and 82 protons as a function of 02- Here, fi\ = 
and 7 = 0. Thomas Fermi energies, second-order corrections and the fourth-order corrections are 
shown in upper, middle and bottom panels respectively. 
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are pulled apart and Coulomb repulsion decreases, thereby making the system more bound. 
The /?4 deformation also affects the proton distribution, but, as expected, the effect of hex- 
adecapole deformation is less prominent in comparison with that of quadrupole deformation. 
Thus, the repulsion among protons does decrease with increasing /3 4 , but the decrease is not 
large enough to make the system more bound with larger f3^. 

By keeping /3 2 and (3± fixed, if the parameter 7 is varied, then it is found that the resulting 
energies are independent of the sign of 7. Moreover, the 7 dependence of the WK energies 
is found to be rather weak. Therefore, here we do not present these result explicitly. 

The fourth-order calculation for protons is very time consuming. Typically, it takes tens of 
minutes to do a complete WK calculation. Most of the run-time being consumed by particle 
number determination and the fourth-order calculations for protons. Thus, it is necessary 
to find an accurate approximation scheme for the fourth-order calculation for protons. Since 
in the nuclear interior, the Coulomb potential has approximately a quadratic nature (see 
Fig. (1)), it is expected that the Coulomb potential will have small influence on the fourth- 
order calculations (note that one needs higher-order derivatives in the fourth-order energy 
calculations). One may therefore drop the Coulomb potential completely from the fourth- 
order corrections; we shall refer to this approximation as "quadratic approximation". This 
approximation has been checked explicitly by performing exact fourth-order calculations 
for protons. The maximum difference between the WK energies obtained by using exact 
calculation and the quadratic approximation is found to be of the order 100 keV for 82 
protons. It turns out that the difference between the quadratic approximation and exact 
calculation decreases with decreasing charge number. This approximation can be improved 
by keeping the Laplacian of the Coulomb potential in the fourth-order contribution i.e., the 
terms of the form (V 2 V) 2 and V 4 V in Eq. ffTTj) . This means that for protons, only the term 
V 2 (VV^) 2 is dropped from Eq. (|T7|) . It is found that with this modification, the value of 
the fourth-order correction energy for the mean field part for protons almost coincides with 
the value obtained by taking all of the derivatives of the Coulomb potential into account. 
This helps in reducing the total runtime further. Thus, effectively, with the interpolation for 
Coulomb potential as discussed before (see section IV), and the approximations introduced 
in the fourth-order correction terms for protons in the present section, the runtime reduces 
from tens of minutes to just about two minutes, without affecting the desired accuracy of 
the calculations. 
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VI. WIGNER-KIRKWOOD SHELL CORRECTIONS AND COMPARISON WITH 
STRUTINSKY CALCULATIONS 



Numerically, it has been demonstrated that the WK and Strutinsky shell corrections 
are close to each other This is expected, since it has recently been shown j|3(3] that 

the Strutinsky level density is an approximation to the semi-classical WK level density 
For illustration, we present and discuss the WK and the corresponding Strutinsky shell 
corrections for the chains of Pb, Gd and Dy isotopes. For the sake of completeness, we first 
present and discuss the essential features of the Strutinsky smoothing scheme. 

According to the Strutinsky smoothing scheme, the smooth level density for a one-body 



Hamiltonian is given by 



49]: 



^^^--'"''^^.(t)' < 44 > 

where q are the single-particle energies calculated by diagonalising the Hamiltonian matrix. 
The smoothing constant 7 is taken to be of the order of huo (hu)o = 1.2 x 41v4~ 1//3 ). N s 
is the smoothing order, and is assumed to be equal t o 6 in the present work; Hj are the 



Hermite polynomials; and Sj is a constant, defined as 

2i (772)! 



49|: 



= O0V ■ /r»M ' 3 eVen ) 



0, for j odd. (45) 



The Strutinsky shell correction is given by: 

Est = J^ei - / eg st {e)de, (46) 
i=i 

where N n is the number of nucleons. This, upon substituting the expression for g st , yields 

H, 



t=i 3=0 



2 L v Jn 2v^F 



72 iV» 



where 



e ^ ^ — \ [7 _ _ 1 I 

—j= Sk \_2 Hk ^ u j) + e i H k~i( u j) + k-fH k _ 2 (uj)^ > , (47) 
v n k=i "~ J 



Uj = ^— ^. (48) 

7 



20 



Here, A is the chemical potential, calculated iteratively from the particle number condition. 
The error integral erf (a;) is defined as: 

erf (a) = 4= / e~ z2 dz. (49) 
V 71 " Jo 

It should be noted that the Strutinsky procedure described here uses positive energy 
states generated by diagonalizing the Hamiltonian matrix, and not by taking resonances 
into account and smoothing them. Further, in practice, the summations defined above do 
not extend up to infinity, but are cut off at a suitable upper limit. The limit is chosen in 
such a way that all the states up to ~ 4hu are included in the sum. It has been shown that 
the uncertainty in the Strutinsky shell corrections obtained this wa y is typically of the order 
of 0.5 MeV [49j. For lighter nuclei, however, it has been concluded 49j that this uncertainty 
is larger. 

The total WK shell correction for the chain of even even Lead isotopes ( 178_214 Pb) is 
plotted in Fig. (TJJ, along with the corresponding values obtained by using the Strutinsky 
smoothing method. It is found that both the WK and Strutinsky results exhibit very 
similar trends. As expected, there is a prominent minimum observed for 208 Pb, indicating 
the occurrence of shell closure. The WK and Strutinsky shell corrections slightly differ from 
each other. The difference is not a constant, and is found to be increasing slowly towards 
the more neutron deficient Lead isotopes. 

Next we plot the calculated (WK) and the corresponding Strutinsky shell corrections 
for the chains of even even Gd and Dy isotopes, with neutron numbers ranging from 72 
to 92. Apart from 144 > 146 . 148 Gd and 146 > 148 > 150 Dy, the rest of the nuclei considered here are 
known to be deformed [l^ . 
the Moller - Nix compilation 
shell corrections agree with each other, within few hundred keVs. The prominent minimum 
at shell closure at neutron number 82 is clearly seen. In these cases as well, the difference 
between the two calculations is not a constant. It is larger in the neutron deficient region, 
and becomes smaller as neutron number increases. 



? or this test run, we adopt the deformation parameters from 
l^ . It is seen that the WK and the corresponding Strutinsky 
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VII. LIPKIN-NOGAMI PAIRING MODEL 



The pairing correlations, important for the open shell nuclei, are often taken into account 
within the framework of the BCS model. The BCS model, however, has two serious short- 
comings: 1) particle number fluctuation (the BCS wavefunctions are not particle number 
eigenstates), and 2) there may exist critical values of the pairing strength, below which the 
BCS equations may not have any non-trivial solutions. In order to overcome these diffi- 
culties, Lipkin, Nogami and co-workers proposed to minimise the expectation value of the 



model Hamiltonian 



29 



30. 



3l|: 



n 



H - AiiV - A 2 iV 2 



(50) 



by determining Ai and A 2 using certain conditions. Here, H is the pairing Hamiltonian, and 
N is the particle number operator. Minimisation of the expectation value of H — \\N leads 
to the usual BCS model, with Ai determined from the particle number condition. Thus, in 
Eq. ([ST)]) above, the quantity Ai is a Lagrange multiplier, but the particle number fluctuation 
constant A 2 is not. 

In practice, the LN calculation is carried out by assuming a constant pairing matrix 
element, G. For a given nucleus (assumed to be even-even for simplicity), one considers 
Nh doubly degenerate states below, and N p doubly degenerate states above the Fermi level. 
These states contain M nucleons. In practice, one takes = N p = N/2 or Z/2, depending 
on whether it is being applied to neutrons or protons. The occupation probabilities v%, the 
pairing gap A, the chemical potential A (= Ai + 2 A? (A/* + 1), see Ref. [3lJ), and the constant 



A 2 are determined iteratively using the conditions 

A = 



13 



3l|: 



such that 



G 22 U k v k 
k 



£k ~ A 



(51) 
(52) 



{(e k -X) 2 -Ai} 



1/2 



(53) 



and 



e k = E k + (4A 2 — G)v 



k ! 



(54) 
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where Ek are the single-particle energies and u\ = 1 — v\. The particle number fluctuation 
constant A2 is given by: 



4 4 

k U k v k 



X = 

The pairing matrix element G is calculated by the Moller-Nix prescription 



(55) 



2 



p L ln|^/ai + A 2 + a 2 | - p L In j^a 2 + A 2 + ai j (56) 



Here, p~£ = gwx/2 is the Wigner-Kirkwood averaged level density (see Eq. ( JTBT) . Factor of 2 
appears because each quantal level here has degeneracy of 2. The level density is evaluated 
at fermi energy); q 2 = N /2p L and a x = —M/2p L and A is the average pairing gap, taken 
to be 3.3/Af 1 / 2 [13]. 



The ground-state energy within the LN model is given by: 

E g = 2^£ fe -^-G£^_ 4A2 £^2 (57) 

fc k k 

The pairing correlation energy, -E pa i r is obtained by subtracting the ground-state energy in 
absence of pairing from Eq. (j57j) : 

= E g -2"£E k - GM/2 . (58) 



VIII. CALCULATION OF BINDING ENERGIES 

As an illustrative example, we now present and discuss the calculated binding energies 
(in this paper, we take binding energies as negative quantities) for 367 even-even, even- 
odd, odd-even and odd-odd spherical nuclei. These nuclei are predicted to be spherical or 
nearly spherical (/5 2 < 0.05) in the Moller-Nix calculations [3] and include 38_52 Ca, 42_54 Ti, 
ioo-i34g n ^ anc j i78-2i4pk rj^g (bailed list of nuclei considered in the present fit can be found 
m Ref. ([51]). Of course, it is known that the prediction of sphericity does depend to some 
extent on the details of the density functional employed [52J. Therefore, it may so happen 
that some of the nuclei assumed to be spherical here, may actually turn out to be slightly 
deformed when energy minimization is carried out on the grid of deformation parameters. 

Our calculation proceeds in the following steps. For each nucleus, the quantum mechani- 
cal and WK energies are calculated as described earlier. This then yields values of the shell 
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corrections (SE) for these nuclei. 



Lipkin-Nogami scheme 
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?he pairing energies (E pair ) are then calculated using the 
31] described previously in the same potential well where 
the shell correction is computed. These two pieces constitute the microscopic part of the 
binding energy. The macroscopic part of the binding energy (El dm) is obtained from the 
liquid drop formula. Thus, for a given nucleus with Z protons and N neutrons (mass number 
A = N + Z), the binding energy in the Mic-Mac picture is given by: 



E(N,Z) = E LDM + 6E + E. 



pair 



(59) 



The liquid drop part of binding energy is chosen to be: 



E 



LDM 



a v 



1 + ^T Z (T Z + 1) 



A + a s 



4/L 



! + ^T Z (T Z + 1) 



A 2/3 



3Z 2 e 2 

+ ■= 7T7T + 



C 4 Z 2 



(60) 



5r AV3 1 A 

where the terms respectively represent: volume energy, surface energy, Coulomb energy 
and correction to Coulomb energy due to surface diffuseness of charge distribution. The 
coefficients a v , a s , k v , k s , r and C4 are free parameters; T z is the third component of 
isospin, and e is the electronic charge. The free parameters are determined by minimising 
the x 2 value in comparison with the experimental energies: 



X 



n 



E 

j=o L 



E{N, V Z 3 )-El 



(J) ■ 

expt 



AE {j) 



(61) 



where E(Nj,Zj) is the calculated total binding energy for the given nucleus, E^J pt is the 
corresponding experimental value 53], and AE^J pt is the uncertainty in E~J t . In the present 
fit, for simplicity, AE^ 3 ] pt is set to 1 MeV. The minimisation is achieved using the well-known 



Levenberg-Marquardt algorithm 54. 155]. 



TABLE I: Values of the liquid drop parameters obtained through the \ 2 minimisation. 



Quantity 


Value 


Quantity 


Value 


a v 


-15.841 (MeV) 


a s 


19.173 (MeV) 


k v 


-1.951 


ks 


-2.577 


ro 


1.187 (fm) 


c 4 


1.247 (MeV) 
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FIG. 5: Difference between the calculated and the experimental [53j binding energies. The corre- 



sponding differences obtained for the Moller-Nix values of binding energies 



131 ] are also presented. 



For the set of nuclei considered here, the rms deviation in binding energies turns out to 
be 630 keV, which, indeed is gratifying. The rms deviation obtained for the same nuclei 
with the Moller-Nix mass formula turns out to be 741 keV. The liquid drop parameters are 
presented in Table I. Clearly, the obtained values of the parameters are reasonable. The 
detailed table containing the nuclei considered in the present fit, and the corresponding 
calculated and experimental 53j binding energies may be found in Ref. 51]. 



To examine the quality of the fit 
and the corresponding experimental 



urther, first, we plot the difference between the fitted 



531 ] binding energies for the 367 nuclei as a function of 



;he mass number A in Fig. (jSj) The corresponding differences obtained for the Moller-Nix 
13] ] values of binding energies are also plotted in the same figure for comparison. It is amply 
clear from the figure that the fitted binding energies are close to the experiment (within 1 
MeV). Overall, the quality of the present fit is slightly better than that of the Moller-Nix fit 
(the rms deviations, respectively, are 630 and 741 keV). Particularly for the lighter nuclei, 
the present calculations are comparatively closer to the experiment. 



Next, the difference between the calculated and the corresponding experimental 



531 ] bind- 



ing energies (denoted by "WK") for Ca, Ti, Sn, and Pb isotopes considered in this fit are 
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presented in Fig. ([6]). The differences obtained by using the Moller-Nix [13( values of bind- 
ing energies (denoted by "MN") are also shown there for comparison. It can be seen that 
the present calculations agree well with the experiment. It is found that the differences vary 
smoothly as a function of mass number: the exceptions being the doubly closed shell nuclei 
48 Ca, 132 Sn, and 208 Pb, where a kink is observed. The overall behaviour of the differences is 
somewhat smoother than that obtained by using the values of Moller and Nix. 
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FIG. 6: Difference between calculated binding energies and experiment [53]. Results are shown 
for the present calculation (WK), for the Moller-Nix values (MN), and using the Rost parameters 
in the Woods-Saxon form factors as described in the text. 



To investigate the effect of the parameters of the single-parti cle p otential, we make a refit 
of the liquid drop parameters, by using the Rost parameters 56( in the microscopic part 
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energies for Sc, Sn and Pb isotopes. 
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of the binding energy computed in the WK approximation. That is, we calculate the shell 
corrections and pairing energies employing the Rost parameters for the Woods-Saxon form 
factors, and then fit the liquid drop parameters again for the same set of nuclei. The rms 
deviation obtained in this case (1.14 MeV for even even nuclei) is much worse that the one 
obtained for the parameters mentioned in Section 3 (the rms deviation obtained for this 
potential is around 0.680 MeV), which is amply clear from the figure. It is well known that 
the Rost parameters have very large half density radii. As a consequence, the values of the 
moment of inertia and rms radii obtained by using these parameters deviate strongly from 
the corresponding experimental values. On the contrary, the parametrisation used in the 
present analysis (see Section 3), yields reasonable values of moments of inertia and radii. 
Thus, overall, this potential is more realistic than the Rost potential. This is reflected in 
the calculated binding energies as well, showing clearly that the choice of the single-particle 
potential (or in other words, the parameters) is indeed important for reliable predictions of 
binding energies (and hence the masses). 

Single and two neutron separation energies (Si n and S2 n ) are crucial observables. They 
are obtained by calculating binding energy differences between pairs of isotopes differing 
by one and two neutron numbers, respectively. The single neutron separation energies 
govern asymptotic behaviour of the neutron density distributions [stJ. They exhibit odd- 
even staggering along an isotopic chain, indicating that the isotopes with even number of 
neutrons are more bound than the neighbouring isotopes with odd number of neutrons. The 
systematics of S2 n primarily reveals the shell structure in an isotopic chain. The correct 
prediction of these separation energies is crucial for determination of the neutron drip lines. 
The calculated Si n and S 2n values for Sc, Sn and Pb isotopes are displayed in Fig. (J7J). The 
corresponding experimental values of Si„ and S^n |53| are also plotted for comparison. The 
agreement between calculations and experiment is found to be excellent. The odd - even 
staggering is nicely reproduced. The shell closures at 132 Sn and 208 Pb are clearly visible 
both in single and two neutron separation energies. At a finer level, however, a marginal 
underestimation of the shell gap at the neutron number 82 (126) is observed in 132 Sn ( 208 Pb). 
Finally, we remark that the calculated single and two proton separation energies are also 
found to be in close agreement with the experiment. 

The results presented in this section indicate that the present calculations of binding 
energies, indeed, are reliable. 
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IX. SUMMARY AND FUTURE OUTLOOK 



In the present work, we intend to carry out reliable mass calculations for the nuclei span- 
ning the entire periodic table. For this purpose, we employ the 'microscopic-macroscopic' 
framework. The microscopic component has two ingredients: the shell correction energy 
and the pairing energy. The pairing energy is calculated by using the well-known Lipkin- 
Nogami scheme. To average out the given one-body Hamiltonian (and hence find the shell 
corrections, given the total quantum mechanical energy of the system), we use the semi- 
classical Wigner-Kirkwood expansion technique. This method does not use the detailed 
single-particle structure, as in the case of the conventional Strutinsky smoothing method. 
In addition to the bound states, the Strutinsky scheme requires the contributions coming 
in from the continuum as well. Treating the continuum is often tricky, and in most of the 
practical calculations, the continuum is taken into account rather artificially, by generating 
positive energy states by means of diagonalisation of the Hamiltonian matrix. For neutron- 
rich and neutron-deficient nuclei, the contribution from the continuum becomes more and 
more important as the Fermi energy becomes smaller (less negative). Uncertainty in the 
conventional Strutinsky scheme thus increases as one goes away from the line of stability. 
It is therefore expected that the Wigner-Kirkwood method will be a valuable and suitable 
option especially for nuclei lying far away from the line of stability. 

We now summarise our observations and future perspectives: 

1. Semi-classical averaging of a realistic one-body Hamiltonian using the Wigner- 
Kirkwood expansion of the partition function correct up to fourth order in h is carried 
out for the deformed systems, both for protons and neutrons. The spin-orbit as well 
as Coulomb potentials are explicitly taken into account. 

2. The smooth energies thus obtained are investigated in detail as a function of three 
deformation parameters: fa, At, and 7. As expected, the energies corresponding to 
the leading-order term in the expansion as well as the correction terms vary smoothly 
as a function of deformation parameters. 

3. Differences between the quantum mechanical and the corresponding averaged energies 
yield the shell corrections. These, along with the pairing energies obtained by using 
the Lipkin-Nogami scheme constitute the "microscopic" part of the nuclear binding 
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energy in the 'Mic-Mac' picture. Using a simple liquid drop ansatz with six adjustable 
parameters, it is demonstrated that the present approach indeed, is feasible, and very- 
promising. For the test case presented here, comprising of 367 spherical nuclei, the 
rms deviation of the predicted binding energies from the experimental values turns 
out to be 630 keV. 

4. The importance of the one-body potential in reliable estimations of nuclear binding 
energies is explicitly demonstrated. It should be noted that the Woods-Saxon param- 
eters used in this work have been fitted for the Coulomb potential calculated by using 
the uniform density (sharp surface) approximation. The Coulomb potential we use 
is obtained from folding Woods-Saxon density profile with the Coulomb interaction. 
Therefore, before performing the large scale calculations, we intend to make a refit to 
the Woods-Saxon potential, with the Coulomb potential obtained from folding. 

5. Having established the feasibility of the present approach, we now intend to extend our 
binding energy calculations to deformed nuclei. For this purpose, we plan to minimise 
the binding energy on a mesh of deformation parameters to find the absolute minimum 
in the deformation space. Work along these lines is in progress. 

APPENDIX A: GEOMETRY OF DISTANCE FUNCTION 

Consider an arbitrary surface, defined by the relation r = r s , where r s is given by Eq. 
( I2TT) of the text. Let us fix the origin of the coordinate system at the centre of mass of the 
object. Let r = (r,9,(j)) define an arbitrary point in space. This point could be inside or 
outside the surface. Here, for concreteness, we assume that it is within the volume of the 
object. Our aim is to find the minimum distance of the point r to the surface r = r s . 
To achieve this, construct a vector R s from the centre of mass to the surface. To find the 
minimum distance, one has to minimise the object \r — R s \. Denoting the angle between the 
r and R s vectors as we have: 



f-R 



■s 



^R 2 S + r 2 - 2rR s cos ^ , 



(Al) 



where, from Fig. (jXJ), the cosine of the angle \1/ is given by: 




(A2) 



31 



FIG. 8: Geometry of distance function 



The latter result can be proved easily by considering a unit sphere, and a spherical triangle 
constructed with unit vectors f, R s , and z. For spherical symmetry, the vectors r and R s 
are parallel to each other, and one recovers the usual spherical Woods-Saxon form factor. 



APPENDIX B: COULOMB POTENTIAL AND ITS DERIVATIVES 
1. Proof of Eq. ([43]) 

The Coulomb potential for an arbitrary charge distribution is given by: 

V c (r) = e 2 [ p^jj—rd? . (Bl) 



Let, for brevity, \r — r*| = 1Z. Consider: 

v> 



K 



Here, the symbol means that the differentiation is done with respect to the r', 9 
coordinates. Let us consider the above derivative component-wise. The contribution coming 
from the first component is: 

d <~R~ ~ K + (B2) 
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Adding contributions coming from all the three components, one gets: 



K 



n 



With this, the potential becomes: 



— e 



K 



d? 



p(f y )V r -v • TZdr 1 



(B3) 

(B4) 
(B5) 



Here, the term in the curly brackets has been represented as a unit vector 1Z. Using the 
identity: 



- Tl , (B6) 
one obtains: 

V c [f) = y J p{7)V 2 F \f-7\d7 (B7) 
which, upon integrating by parts and transferring derivatives to density, becomes: 

2 f 

V c (r) = e — j df*\r - f|V£p(0 (B8) 

q.e.d. 



2. Derivatives of Coulomb Potential 

The calculation of the higher-order derivatives of the Coulomb potential (third and 
above), even with the form defined in Eq. (1431) . turns out to be numerically unstable. 
For this purpose, we employ the Poisson's equation. According to this, the Laplacian of 
Coulomb potential is proportional to the charge density: 

VV c (f) = -4vre 2 p(r) (B9) 

The Laplacian of V 2 Vc{r) is simple to compute, for, all one needs to calculate there are the 
derivatives of density (assumed to be of Woods-Saxon form). 

Thus, it is desirable to generate the required higher order-derivatives of the Coulomb 
potential (see expression (T42l) in the text) from Poisson's equation. For this purpose, we 
evaluate the commutators: 

' 2 ld Id , 
V , V 

r 86 r 86 



„J 8 „ 
V — , — V 

8r 8r 



Vb(r), 



Vb(r), 



9 esc 6 8 esc 6 8 



r 86' r 86 



Vc(r) 
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The results are: 



9 d 
or 



1 d 



d 9 2 
F W C + - 
or r 



' 13 a 2 



■7- l'^|y C 
r otp 



1 <9 _ 9tj . esc 2 
csc# a 



§e Vc + 2cote W Vc 



2 a 2 



90 



„, T , csc 3 # <9 T , 2csc# 
V 2 Vc + — 3-^-rVc =- 



r 2 dOdr 
cotO d 2 



V c 



(BIO) 
(Bll) 



d9d(f) 



d 2 ' 

Vc + -^tt-Vc 



dOdr 



(B12) 



With these expressions, the required higher-order derivatives of the Coulomb potential can 
be generated. These are then used to evaluate the fourth-order WK energy, as we have 
described in Section 4. 



ACKNOWLEDGMENTS 

AB is thankful to KTH Stockholm and IPN Orsay for financial support. M.C. and 
X.V. partially supported by the Consolider Ingenio 2010 Programme CPAN CSD2007-00042 
and grants FIS2008-01661 from MEC and FEDER and 2009SGR-1289 from Generalitat de 
Catalunya. 



[1] D. Seweryniak and T. L. Khoo (eds.), Nuclei at the Limits, Argonne, USA, 26-30 July 2004, 
AIP Conference Proceedings, Vol. 764 (Springer, 2005). 

[2] Proceedings of the Conference ENAM'04 on Exotic Nuclei and Atomic Masses, Pine Mountain, 
USA, 12-16 September 2004, edited by C. J. Gross, W. Nazarewicz, and K. P. Rykaczewski, 
Eur. Phys. J. A 25, sOl (2005); Proceedings of the Conference ENAM'08 on Exotic Nuclei and 
Atomic Masses, Ryn, Poland, 7-13 September 2008, to be published in Eur. Phys. J. A. 

[3] D. Lunney, J. M. Pearson and C. Thibault, Rev. Mod. Phys. 75, 1021 (2003). 

[4] K. Blaum, Phys. Rep. 425, 1 (2006). 

[5] Yu.Ts. Oganessian, JINR Reprint E7-2002-64. 

[6] Yu.Ts. Oganessian et al, Nucl. Instr. and Meth. B204, 606 (2003). 

[7] V. S. Kolhinen et al, Nucl. Instr. and Meth. A528, 776 (2004). 

[8] S. Goriely et al, Phys. Rev. C 66, 024326 (2002) and references cited therein. 



34 



[9] S. Goriely, N. Chamel and J. M. Pearson, Phys. Rev. Lett. 102, 152503 (2009). 
[10] S. Hilaire and M. Girod, International Conference on Nuclear Data for Science and Technology, 



[11 
[12 

[13 
[14 
[15 
[16 
[17 
[18 
[19 
[20 
[21 
[22 
[23 
[24 
[25 
[26 

[27; 

[28 
[29 
[30 
[31 
[32 
[33 
[34 
[35 
[36 



2007 (DOI:10.1051/ndata:07709). Also, see |http://www-phynu.ce^ /HFB-Gogn yhtni] for the 
resulting 'AMEDEE' structure database. 

G. A. Lalazissis, S. Raman and P. Ring, At. Data Nucl. Data Tables 71, 1 (1999). 

P. Moller, J. R. Nix, W. D. Myers and W. J. Swiatecki, At. Data Nucl. Data Tables 59, 185 

(1995). 

P. Moller, J. R. Nix and K. -L. Kratz, At. Data Nucl. Data Tables 66, 131 (1997). 

K. Pomorski and J. Dudek, Phys. Rev. C 67, 044316 (2003) and references cited therein. 

W. D. Myers and W. J. Swiatecki, Nucl. Phys. A601, 141 (1996). 

V. M. Strutinsky, Nucl. Phys. A95, 420 (1967). 

G. G. Bunatian, V. M. Kolomietz and V. M. Strutinsky, Nucl. Phys. A188, 225 (1972). 

E. Wigner, Phys. Rev. 40, 749 (1932). 

J. G. Kirkwood, Phys. Rev. 44, 33 (1933). 

B. K. Jennings, R. K. Bhaduri and M. Brack, Nucl. Phys. A 253, 29 (1975). 

P. Ring and P. Schuck, The Nuclear Many-Body Problem (Springer- Ver lag, Berlin, 1980). 

M. Brack, C. Guet and H. -B. Hakansson, Phys. Rep. 123, 275 (1985). 

M. Brack and R. K. Bhaduri, Semi-classical Physics (Addison - Wesley Publishing Co., 1997). 
M. Centelles et al, Phys. Rev. C 74, 034332 (2006) and references cited therein. 
M. Centelles, P. Schuck and X. Vihas, Ann. Phys. 322, 363 (2007). 
S. Shlomo, Nucl. Phys. A539, 17 (1991). 

W. Nazarewicz, T. R. Werner and J. Dobaczewski, Phys. Rev. C 50, 2860 (1994). 
T. Vertse et al, Phys. Rev. C 57, 3089 (1998). 

H. J. Lipkin, Ann. Phys. (N.Y.) 9, 272 (1960). 
Y. Nogami, Phys. Rev. 134, B313 (1964). 

H. C. Pradhan, Y. Nogami and J. Law, Nucl. Phys. A201, 357 (1973). 

F. R. Xu, R. Wyss, and P.M. Walker, Phys. Rev. C 60, 051301 (R) (1999). 

B. K. Jennings, R. K. Bhaduri and M. Brack, Phys. Rev. Lett. 34, 228 (1975). 

PSchuck and X. Vihas, Phys. Lett. B302, 1 (1993). 

M. Centelles et al, Ann. Phys. 266, 207 (1998). 

B. Grammaticos and A. Voros, Ann. Phys. 123, 359 (1979). 

35 



[37; 

[38 
[39 
[40 
[41 
[42' 
[43' 
[44 
[45 
[46 
[47; 
[48 
[49 
[50 
[51 



B. Grammaticos and A. Voros, Ann. Phys. 129, 153 (1980). 

M. Centelles et al, Nucl. Phys. A510, 397 (1990). 

R. D. Woods and D. S. Saxon, Phys. Rev. 95, 577 (1954). 

J. Dudek et al, J. Phys. G: Nucl. Phys. 5, 1359 (1979) 

J. Dudek, Z. Szymahski and T. Werner, Phys. Rev. C 23, 920 (1981) 

W. Nazarewicz et al, Nucl. Phys. A435, 397 (1985). 

S. Cwiok et al, Comp. Phys. Comm. 46, 379 (1987). 

R. Wyss et al, Nucl. Phys. A511, 324 (1991). 

J. Damgaard et al, Nucl. Phys. A135, 432 (1969). 

J. Bennewicz and P. K. Huag, Z. Phys. 212, 212 (1968). 

R. Wyss, unpublished. 

W. Satula and R. Wyss, Rep. Prog. Phys. 68, 131 (2005) and references cited therein. 
M. Bolsterli et al, Phys. Rev. C 5, 1050 (1972) and references cited therein. 
B. Mohammed-Azizi and D. E. Medjadi, Phys. Rev. C 74, 054302 (2006). 
For the table, see EPAPS document associated to this article, 



URL http ://www.aip.org/pubservs/epsps.htmH 



The table can also be found at http://www.ecm.ub.es/~xavier/be_WK.dat 



Notice that we take binding energies as negative quantities. 



[52] J. Terasaki, J. Engel and G. F. Bertsch, [irXivi0801.2346v2 [nucl-th]. 
[53] A. H. Wapstra, G. Audi and C. Thibault, Nucl. Phys. A 729, 129 (2003). 
[54] D. W. Marquardt, J. Soc. Ind. Appl. Math. 11, 431 (1963). 

[55] W. H. Press, S. A. Teukolsky, W. H. Vetterling and B. P. Flannery, Numerical Receipes in 

FORTRAN (2 nd edition) (Cambridge University Press, 1992), p. 678. 
[56] E. Rost, Phys. Lett. B26, 184 (1968). 

[57] Y. K. Gambhir and S. H. Patil, Z. Phys. A321, 161 (1985); A. Bhagwat, Y. K. Gambhir and 
S. H. Patil, Eur. Phys. J. A8, 511 (2000). 



36 



